#-*- coding=utf-8 -*-
import partie1 as p1
import numpy as np
import matplotlib.pyplot as plt

def solve_first_eq(b, nmax, method = p1.step_kutta):
    "Résout la première équation pour une valeur donnée de b"
    t0 = np.array([[0.]])
    y0 = np.array([[1.],[0.], [0.]])
    a = 10.
    c = 200.
    d = 50.
    tmax = 1.
    f = lambda X,T: np.array([[a + b*X[2][0] - X[0][0] - X[0][0]*X[1][0]*X[1][0]],
                              [c*(X[0][0] + X[0][0]*X[1][0]*X[1][0] - X[1][0])],
                              [d*(X[1][0] - X[2][0])]])
    (res, RES) = p1.meth_n_step (y0, t0, nmax, tmax/nmax, f, method)
    Xaxis = np.arange(0.,tmax*(1.+1./nmax),tmax/nmax)
    (Xaxis,res) = periodic_value(Xaxis,res)
    x = []
    y = []
    z = []
    for i in range(len(res)):
        x.append(res[i][0][0])
        y.append(res[i][1][0])
        z.append(res[i][2][0])
    g = plt.plot(Xaxis, x)
    h = plt.plot(Xaxis, y)
    i = plt.plot(Xaxis, z)
    plt.legend((g,h,i),("X(t)","Y(t)","Z(t)"))
    plt.xlabel("Temps")
    plt.ylabel("Concentration")
    plt.show()

def periodic_value(x,res):
    "Retourne la partie périodique du tableau res"
    x0 = []
    res0 = []
    m = 0.
    for i in range(len(x)):
        if(x[i] > 0.3 and abs(res[i][0][0] + res[i][0][0]*res[i][1][0]*res[i][1][0] - res[i][1][0])<0.01):
            m = res[i][1][0]
            j = i+20
            while((abs(res[j][1][0]-m)>0.1 or abs(res[j][0][0] + res[j][0][0]*res[j][1][0]*res[j][1][0] - res[j][1][0])>0.01) and j<len(res)-1):
                x0.append(x[j])
                res0.append(res[j])
                j = j+1
            return (x0,res0)
    return (x0,res0)

def val_max(res):
    "Retourne le nombre de maxima locaux de Y"
    val = []
    for i in np.arange(1,len(res)-1,1):
        if(res[i-1][1][0] < res[i][1][0] and res[i][1][0] > res[i+1][1][0]):
            val.append(res[i][1][0])
    return val

def bifurcation(method = p1.step_kutta):
    "Trace le diagramme de bifuraction du système"
    t0 = np.array([[0.]])
    y0 = np.array([[1.],[0.], [0.]])
    a = 10.
    c = 200.
    d = 50.
    tmax = 0.7
    nmax = 7000
    for b in np.arange(0.1,0.2,0.005):
        print b*1000.-100.,"%"
        f = lambda X,T: np.array([[a + b*X[2][0] - X[0][0] - X[0][0]*X[1][0]*X[1][0]],
                              [c*(X[0][0] + X[0][0]*X[1][0]*X[1][0] - X[1][0])],
                              [d*(X[1][0] - X[2][0])]])
        res = p1.meth_n_step (y0, t0, nmax, tmax/nmax, f, method)[0]
        Xaxis = np.arange(0.,tmax*(1.+1./nmax),tmax/nmax)
        res = periodic_value(Xaxis,res)[1]
        val = val_max(res)
        if(len(val) == 3):
            print b
        for i in range(len(val)):
            plt.plot(b,val[i],',',color="black")

solve_first_eq(0.155, 10000)
print "Début du tracé du diagramme..."
bifurcation()
print "Tracé du diagramme terminé."
